---
title: "ReshoringReplicationFinal"
author: "Ricky Clark, Adi Rao"
date: "2024-12-20"
output: html_document
---


```{r Preliminaries}
### Preliminaries
setwd("/Users/USER/FOLDER/SUBFOLDER/Reshoring")
set.seed(7777777)
library(pacman)
pacman::p_load(dplyr, ggplot2, quanteda, plotly, foreign, interplot, plm, reshape2, lda,
               countrycode, sandwich, lmtest, MASS, RColorBrewer, states, mice, VIM, mediation, margins, clusterSEs, optimx, coefplot, systemfit, randomizr, estimatr, ri2, stargazer, fastDummies, cowplot, egg, cobalt, plotrix, tm, stm, stargazer, Rmisc, modelsummary, gridExtra, readr, readxl, broom)

# Wave 1 data
test <- read.csv("reshoring.csv")

# Wave 2 data
results <- read.csv("reshoring_qualtrics_num.csv")

## The globalization variable was coded incorrectly in Wave 2: survey answers 3 and 4 were blended together. We recoded it in Qualtrics, redownloaded it, and appended the corrected variable into the results df. (Does not change data integrity, only fixes leveling of the globalization variable.)
glob_only <- read.csv("reshoring_qualtrics_num2.csv")

glob_only <- glob_only %>% select(globalization)

glob_only <- rename(glob_only,
                       correctglob="globalization")

results <- cbind(results, glob_only)


# For semiconductor chips in the news graph
news <- read.csv("ChipsNews.csv")



```

``` {r Wave 1 Analysis Code}
#############################
#### Survey Wave 1 ##########
#############################

#### Key variable coding ####
test <- test[test$consent_check == 1,] # Delete those that fail consent
test$china <- ifelse(!is.na(test$CHSemi) | !is.na(test$CHSteel), 1, 0) # code China treatment
test$southkorea <- ifelse(!is.na(test$SKSemi) | !is.na(test$SKSteel), 1, 0) # code Korea treatment
test$steel <- ifelse(!is.na(test$SKSteel) | !is.na(test$CHSteel), 1, 0) # code steel treatment
test$semi <- ifelse(!is.na(test$SKSemi) | !is.na(test$CHSemi), 1, 0) # code semi treatment

#### Table with N Sizes ####
# Values are calculated below and manually inputted into a LaTeX table
sum(test$semi == 1 & test$china == 1, na.rm=T)
sum(test$steel == 1 & test$china == 1, na.rm=T)
sum(test$semi == 1 & test$southkorea == 1, na.rm=T)
sum(test$steel == 1 & test$southkorea == 1, na.rm=T)

mean(test$dv[test$semi == 1 & test$china == 1], na.rm=T) # 3.86
mean(test$dv[test$steel == 1 & test$china == 1], na.rm=T) # 3.93
mean(test$dv[test$semi == 1 & test$southkorea == 1], na.rm=T) # 3.79
mean(test$dv[test$steel == 1 & test$southkorea == 1], na.rm=T) # 3.92

#### Data Pre-processing ####
test <- test[test$china == 1 | test$southkorea == 1,] # delete incompletes
test[3:43] <- sapply(test[3:43], as.numeric)
test[48:89] <- sapply(test[48:89], as.numeric) # code as numeric
test$open[test$open %in% c("Na", "N/A", "Yes", "No", "1", "2", "3", "4", "5", "81652", "N/a", "Djdjdjd", "n/a",
          "na" ,"Na" ,"yes" ,"no" ,"none", "None", "Nothing", "Ns", "nothing", "Ya yo no yo ya ya I", "Yeah", ".",
          "Bbbbbbb cndlsprotkt vkfoeiwwi fuchsia ur hhohotp hdjwofjtd cfidkdir fit yo rich ghehruriricu hchdjdoeododo h Shawn half hall Harold Camilla sent healing Hann",
          "I am", "Ligma sugma ", "No ", "k", "Nada", "N/a ", "Ion")] <- NA
test$open_qual <- ifelse(is.na(test$open), 0, 1)
# test <- test[test$Duration..in.seconds. > 300,]
test <- test[test$open_qual == 1,]
test$pol_knwldgeagg <-  test$pol_knowledge_1 + test$pol_knowledge_2 # create political knowledge index
test$coop_interagg <- as.numeric(as.character(test$coop_inter_1)) + 
  as.numeric(as.character(test$coop_inter_2)) + 
  as.numeric(as.character(test$coop_inter_3)) + 
  as.numeric(as.character(test$coop_inter_4)) # create cooperative internationalism index
test$coop_interagg <- test$coop_interagg / 4 # convert index to five-point scale
test$isoagg <- as.numeric(as.character(test$iso_1)) + 
  as.numeric(as.character(test$iso_2)) + 
  as.numeric(as.character(test$iso_3)) # create aggregate isolationism measure
test$isoagg <- test$isoagg / 3 # convert index to five-point scale
test$natsupagg <- as.numeric(as.character(test$nat_sup_1)) + as.numeric(as.character(test$nat_sup_2)) # create national superiority index
test$natsupagg <- test$natsupagg / 2 # convert index to five-point scale
colnames(test)[32] <- "age1"
test <- test %>% 
  #rowwise will make sure the sum operation will occur on each row
  rowwise() %>% 
  #then a simple sum(..., na.rm=TRUE) is enough to result in what you need
  mutate(pat = sum(ch_pat_2, st_pat_2, na.rm=TRUE))
test <- test %>% 
  #rowwise will make sure the sum operation will occur on each row
  rowwise() %>% 
  #then a simple sum(..., na.rm=TRUE) is enough to result in what you need
  mutate(values = sum(values_steel_1, values_chip_2, na.rm=TRUE))
test <- test %>% 
  #rowwise will make sure the sum operation will occur on each row
  rowwise() %>% 
  #then a simple sum(..., na.rm=TRUE) is enough to result in what you need
  mutate(pricea = sum(st_price_auto, ch_price_auto, na.rm=TRUE)) %>%
  mutate(pricel = sum(st_price_lap, ch_price_lap, na.rm=TRUE))
test <- test %>% 
  #rowwise will make sure the sum operation will occur on each row
  rowwise() %>% 
  #then a simple sum(..., na.rm=TRUE) is enough to result in what you need
  mutate(accessa = sum(st_access_auto, ch_access_auto, na.rm=TRUE)) %>%
  mutate(accessi = sum(st_access_inf, ch_access_inf, na.rm=TRUE)) %>%
  mutate(accessd = sum(st_access_sec, ch_access_sec, na.rm=TRUE)) %>%
  mutate(friend = sum(SK_friend, CH_friend, na.rm=TRUE))

#### Bootstrapped Treatment Effects ####
# Write function for calculating bootstrapped quantities of interest
bootTreat2 <- function(dat, dv="dv", B=1500){
  bootResults <- matrix(NA, nrow=B, ncol=5)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    dep.var <- eval(parse(text=paste("temp",dv,sep="$")))
    A <- mean(dep.var[which(temp$china==1)], na.rm=TRUE)
    B <- mean(dep.var[which(temp$southkorea==1)], na.rm=TRUE)
    C <- mean(dep.var[which(temp$semi==1)], na.rm=TRUE)
    D <- mean(dep.var[which(temp$steel==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-B)
    bootResults[i,2] <- (C-D)
    drop(list())
  }
  return(list(dv=dv, boot=bootResults, Ch.Effect=mean(bootResults[,1], na.rm=TRUE), Semi.Effect=mean(bootResults[,2], na.rm=TRUE), na.rm=TRUE))
}
set.seed(26)
## ATE with funding DV
#-># Aggregate
prep.full <- bootTreat2(test, dv="dv")
prep.full$Ch.Effect # 0.042
prep.full$Semi.Effect # -0.115
#To calculate p-values:
1 - length(which(prep.full$boot[,1] > 0))/nrow(prep.full$boot) # p = 0.163
1 - length(which(prep.full$boot[,2] < 0))/nrow(prep.full$boot) # p = 0.002

t.test(test$dv[test$china == 1], test$dv[test$southkorea==1]) # p = 0.326
t.test(test$dv[test$semi == 1], test$dv[test$steel==1]) # p = 0.006

#### Descriptive Statistics ####
desc <- test[c(89, 101, 104:107, 91, 94, 99, 24, 80, 33, 35, 27, 74)] # duplicate data
desc <- data.frame(desc)
desc$republican <- ifelse(desc$partisan_dri == 2, 1, 0)
desc$democrat <- ifelse(desc$partisan_dri == 1, 1, 0)
desc$independent <- ifelse(desc$partisan_dri == 3, 1, 0)
desc <- desc[-c(14)]
covs_desc <- c("Reshoring", "American values", "Access automobiles", "Access infrastructure", "Security concerns", 
               "Country thermometer", "China", "Semiconductor", "National superiority", "Male",
               "Age", "Income", "Education", "Trump voter", "Republican", "Democrat", "Independent") # covariate labels
stargazer(desc, omit.summary.stat = c("p25", "p75"), digits = 2, style = "ajps",
          covariate.labels = covs_desc) # print table; presented as Table A3

plot(density(na.omit(test$dv)), main = "Distribution of DV")



#### Regression Modeling ####
#### Main regressions 
reg1 <- lm(dv ~ china + semi, data = test)
reg1a <- lm(dv ~ china + semi + natsupagg + gender_male + age + 
              income + education + I(partisan_dri == 2) + trump, data = test)
stargazer(reg1, reg1a, omit.stat = c("rsq", "ser", "f"),
          covariate.labels = c("China", "Semiconductor", "National superiority",
                               "Male", "Age", "Income", "Education", "Republican",
                               "Trump voter"))
# Results presented in Table 1

#### Support by Treatment ####
B <- 1500 #Number of bootstraps
dens <- data.frame(y=c(prep.full$boot[,1], prep.full$boot[,2]),
                   dv=rep(c("Support for Reshoring"), each=2*B), 
                   Treatment=rep(c("China", "Semiconductor"), each=B))  

m.cost <- mean(dens$y[dens$Treatment == "China"])
se.cost <- sd(dens$y[dens$Treatment == "China"])
ci.up.cost <- m.cost + (se.cost*1.645) 
ci.low.cost <- m.cost  - (se.cost*1.645)

m.inf <- mean(dens$y[dens$Treatment == "Semiconductor"])
se.inf <- sd(dens$y[dens$Treatment == "Semiconductor"])
ci.up.inf <- m.inf + (se.inf*1.645)  
ci.low.inf <- m.inf  - (se.inf*1.645) 

means.fund <- c(m.cost, m.inf)
fund.ci.up <- c(ci.up.cost, ci.up.inf)
fund.ci.low <- c(ci.low.cost, ci.low.inf)
names <- c("China", "Semiconductor")
ciplot <- data.frame(cbind(means.fund, fund.ci.up, fund.ci.low, names))
ciplot$means.fund <- as.numeric(ciplot$means.fund)
ciplot$fund.ci.low <- as.numeric(ciplot$fund.ci.low)
ciplot$fund.ci.up <- as.numeric(ciplot$fund.ci.up)

ggplot(ciplot, aes(y = means.fund, x = names)) +
  theme_bw() + 
  geom_point(size = 6) +
  geom_errorbar(aes(ymax = fund.ci.up, ymin = fund.ci.low),
                width=.2,
                position=position_dodge(.9)) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_text(size = 28), axis.title.y = element_text(size = 28),
        legend.position = "none", axis.text.x = element_text(size = 26),
        axis.text.y = element_text(size = 26)) +
  geom_hline(yintercept=0, linetype='dotted', col = 'red') +
  labs(y = "Point Estimate", x = "Treatment") +
  scale_fill_grey()





### Reviewers advised we should present distribution plots of the DV not as density (which is smooth given integer data). 
test$dv


# Create a temporary factor variable for plotting purposes
test <- test %>%
  mutate(dv_label = factor(dv, levels = 1:5, labels = c("Strongly Disagree", "Disagree", "Neutral", "Agree", "Strongly Agree")))

# Calculate frequencies and percentages
total_count <- sum(table(test$dv_label))
freq_data <- test %>%
  count(dv_label) %>%
  mutate(percentage = n / total_count * 100)

# Define custom colors
custom_colors <- c("Strongly Disagree" = "#636363", "Disagree" = "#969696", "Neutral" = "#bdbdbd", "Agree" = "#d9d9d9", "Strongly Agree" = "#f0f0f0")

# Create bar chart with custom colors
ggplot(freq_data, aes(x = dv_label, y = percentage, fill = dv_label)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = sprintf("%.1f%%", percentage)), vjust = -0.5) +
  scale_x_discrete(name = "Support for Reshoring Policies") +
  scale_y_continuous(name = "Percentage") +
  scale_fill_manual(values = custom_colors, name = "Support Level") + # Custom subdued colors
  theme_minimal() +
  theme(legend.position = "none")
ggsave("/Users/USER/FOLDER/SUBFOLDER/Reshoring/dvdist3.pdf")













#### Support by treatment figures ####
testp <- test
testp$supp <- ifelse(testp$dv == 4 | testp$dv == 5, 1, 0)

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- plyr::rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


plot1 <- summarySE(testp, measurevar="supp", groupvars=c("china"))
plot2 <- summarySE(testp, measurevar="supp", groupvars=c("semi"))
plot1 <- plot1[c(3,5)]
colnames(plot1)[1] <- "Support"
plot1$Treatment <- c("South Korea", "China")
plot2 <- plot2[c(3,5)]
colnames(plot2)[1] <- "Support"
plot2$Treatment <- c("Steel", "Semiconductors")
testplot <- rbind(plot1, plot2)
testplot <- testplot[c(2, 1, 3, 4),]
testplot$num <- round(testplot$Support, 2)
testplot$Treatment <- factor(testplot$Treatment, levels = testplot$Treatment)

ggplot(testplot, aes(y = Support, x = Treatment, fill = c("blue", "blue", "red", "red"))) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_errorbar(aes(ymin=Support-se, ymax=Support+se), width=0.3) +
  geom_text(aes(label=num), position=position_dodge(width=0.9), vjust=-2.5, size = 8) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_text(size = 32),
        legend.position = "none", axis.text.x = element_text(size = 30, angle = 30, vjust = 0.6),
        axis.text.y = element_text(size = 30)) +
  coord_cartesian(ylim=c(0.5,0.8))

plot3 <- summarySE(testp, measurevar="supp", groupvars=c("partisan_dri"))
plot3 <- plot3[c(1,3,5)]
colnames(plot3) <- c("Party", "Support", "se")
plot3$Party <- c("Democrat", "Republican", "Independent")
plot3$num <- round(plot3$Support, 2)
plot3$Party <- factor(plot3$Party, levels = plot3$Party)

ggplot(plot3, aes(y = Support, x = Party, fill = c("red", "blue", "green"))) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_errorbar(aes(ymin=Support-se, ymax=Support+se), width=0.3) +
  geom_text(aes(label=num), position=position_dodge(width=0.9), vjust=-2.5, size = 8) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_text(size = 32),
        legend.position = "none", axis.text.x = element_text(size = 30, angle = 30, vjust = 0.6),
        axis.text.y = element_text(size = 30)) +
  coord_cartesian(ylim=c(0.5,0.8))

#### Support by treatment table ####
mean(test$dv[test$china == 1 & test$semi == 1]) # 3.87
mean(test$dv[test$china == 1 & test$semi == 0]) # 3.95
mean(test$dv[test$southkorea == 1 & test$semi == 1]) # 3.79
mean(test$dv[test$southkorea == 1 & test$semi == 0]) # 3.94


mean(test$dv[test$semi == 1 & test$aware_semi==1])
mean(test$dv[test$semi == 1 & test$aware_semi==0])

#### Availability concerns ####
t.test(test$ch_access_auto,test$st_access_auto) # 3.63 vs. 3.5, p = 0.002
t.test(test$ch_price_auto,test$st_price_auto) # no price difference




```

```{r Wave 2 Analysis Code}

#############################
#### Survey Wave 2 ##########
#############################

#### Table with N Sizes ####
# Values are calculated below and manually inputted into a LaTeX table
sum(results$price==1250 & results$country=="Canada" & results$jobs=="blue collar", na.rm=T) # 312
sum(results$price==1100 & results$country=="Canada" & results$jobs=="blue collar", na.rm=T) # 270
sum(results$price==1250 & results$country=="China" & results$jobs=="blue collar", na.rm=T) # 270
sum(results$price==1100 & results$country=="China" & results$jobs=="blue collar", na.rm=T) # 290
sum(results$price==1250 & results$country=="Canada" & results$jobs=="white collar", na.rm=T) #297
sum(results$price==1100 & results$country=="Canada" & results$jobs=="white collar", na.rm=T) # 282
sum(results$price==1250 & results$country=="China" & results$jobs=="white collar", na.rm=T) # 266
sum(results$price==1100 & results$country=="China" & results$jobs=="white collar", na.rm=T) # 267

mean(results$reshoring[results$price==1250 & results$country=="Canada" & results$jobs=="blue collar"], na.rm=T) # 3.64
mean(results$reshoring[results$price==1100 & results$country=="Canada" & results$jobs=="blue collar"], na.rm=T) # 3.77
mean(results$reshoring[results$price==1250 & results$country=="China" & results$jobs=="blue collar"], na.rm=T) # 3.72
mean(results$reshoring[results$price==1100 & results$country=="China" & results$jobs=="blue collar"], na.rm=T) # 3.92
mean(results$reshoring[results$price==1250 & results$country=="Canada" & results$jobs=="white collar"], na.rm=T) # 3.59
mean(results$reshoring[results$price==1100 & results$country=="Canada" & results$jobs=="white collar"], na.rm=T) # 3.76
mean(results$reshoring[results$price==1250 & results$country=="China" & results$jobs=="white collar"], na.rm=T) # 3.70
mean(results$reshoring[results$price==1100 & results$country=="China" & results$jobs=="white collar"], na.rm=T) # 3.82

#### Descriptive Statistics ####
desc <- results[c(36, 22, 31, 32, 34, 26, 42)] # duplicate data
desc$partisan_dri <- ifelse(desc$partisan_dri == 2, 1, 0)
covs_desc <- c("Reshoring", "Male", "Age", "Income", "Education", "Republican", "Trump voter") # covariate labels
stargazer(desc, omit.summary.stat = c("p25", "p75"), digits = 2, style = "ajps",
          covariate.labels = covs_desc) # print table; presented as Table A3

plot(density(na.omit(results$reshoring)), main = "Distribution of DV")
#### Bootstrapped Treatment Effects ####
bootTreat2 <- function(dat, dv="dv", B=1500){
  bootResults <- matrix(NA, nrow=B, ncol=5)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    dep.var <- eval(parse(text=paste("temp",dv,sep="$")))
    A <- mean(dep.var[which(temp$country=="China")], na.rm=TRUE)
    B <- mean(dep.var[which(temp$country=="Canada")], na.rm=TRUE)
    C <- mean(dep.var[which(temp$jobs=="blue collar")], na.rm=TRUE)
    D <- mean(dep.var[which(temp$jobs=="white collar")], na.rm=TRUE)
    E <- mean(dep.var[which(temp$price==1100)], na.rm=TRUE)
    G <- mean(dep.var[which(temp$price==1250)], na.rm=TRUE)
    bootResults[i,1] <- (A-B)
    bootResults[i,2] <- (C-D)
    bootResults[i,3] <- (E-G)
    drop(list())
  }
  return(list(dv=dv, boot=bootResults, Ch.Effect=mean(bootResults[,1], na.rm=TRUE), 
              Jobs.Effect=mean(bootResults[,2], na.rm=TRUE),
              Price.Effect=mean(bootResults[,3], na.rm=TRUE)))
}
set.seed(26)
## ATE with funding DV
#-># Aggregate
prep.full <- bootTreat2(results, dv="reshoring")
prep.full$Ch.Effect # 0.104
prep.full$Jobs.Effect # 0.046
prep.full$Price.Effect # 0.159
#To calculate p-values:
1 - length(which(prep.full$boot[,1] > 0))/nrow(prep.full$boot) # p = 0.008
1 - length(which(prep.full$boot[,2] > 0))/nrow(prep.full$boot) # p = 0.144
1 - length(which(prep.full$boot[,3] > 0))/nrow(prep.full$boot) # p = 0.0007

t.test(results$reshoring[results$country=="China"], 
       results$reshoring[results$country=="Canada"]) # p = 0.012
t.test(results$reshoring[results$jobs=="blue collar"], 
       results$reshoring[results$jobs=="white collar"]) # p = 0.258
t.test(results$reshoring[results$price==1100], 
       results$reshoring[results$price==1250]) # p = 0.0001


#### Support by treatment figures ####
testp <- results
testp$reshoring <- ifelse(testp$reshoring == 4 | testp$reshoring == 5, 1, 0)
plot1 <- summarySE(testp, measurevar="reshoring", groupvars=c("country"))
plot1 <- plot1[c(1,3,5)]
colnames(plot1) <- c("Treatment", "Support", "se")
plot2 <- summarySE(testp, measurevar="reshoring", groupvars=c("jobs"))
plot2 <- plot2[c(1,3,5)]
colnames(plot2) <- c("Treatment", "Support", "se")
plot3 <- summarySE(testp, measurevar="reshoring", groupvars=c("price"))
plot3 <- plot3[c(1,3,5)]
colnames(plot3) <- c("Treatment", "Support", "se")
testplot <- rbind(plot1, plot2, plot3)
testplot$Treatment[testplot$Treatment=="white collar"] <- "White collar"
testplot$Treatment[testplot$Treatment=="blue collar"] <- "Blue collar"
testplot$num <- round(testplot$Support, 2)
testplot$Treatment <- factor(testplot$Treatment, levels = testplot$Treatment)

ggplot(testplot, aes(y = Support, x = Treatment, fill = c("red", "red", 
                                                          "blue", "blue", 
                                                          "green", "green"))) +
  geom_bar(position = 'dodge', stat = "identity") +
  ylim(0.0, 0.8) +
  geom_errorbar(aes(ymin=Support-se, ymax=Support+se), width=0.3) +
  geom_text(aes(label=num), position=position_dodge(width=0.9), vjust=-2.5, size = 8) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_text(size = 32),
        legend.position = "none", axis.text.x = element_text(size = 30, angle = 30, vjust = 0.6),
        axis.text.y = element_text(size = 30)) +
  coord_cartesian(ylim=c(0.5,0.8))

plot3 <- summarySE(testp, measurevar="reshoring", groupvars=c("partisan_dri"))
plot3 <- plot3[c(1,3,5)]
colnames(plot3) <- c("Party", "Support", "se")
plot3$Party <- c("Democrat", "Republican", "Independent")
plot3$num <- round(plot3$Support, 2)
plot3$Party <- factor(plot3$Party, levels = plot3$Party)


ggplot(plot3, aes(y = Support, x = Party, fill = c("red", "blue", "green"))) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_errorbar(aes(ymin=Support-se, ymax=Support+se), width=0.3) +
  geom_text(aes(label=num), position=position_dodge(width=0.9), vjust=-2.5, size = 8) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_text(size = 32),
        legend.position = "none", axis.text.x = element_text(size = 30, angle = 30, vjust = 0.6),
        axis.text.y = element_text(size = 30)) +
  coord_cartesian(ylim=c(0.5,0.8))

# supp or not
testp <- results
testp$supp <- ifelse(testp$reshoring == 4 | testp$reshoring == 5, 1, 0)
testaggc <- aggregate(testp$supp, list(testp$country), FUN = mean)
colnames(testaggc) <- c("Treatment", "Support")
testaggs <- aggregate(testp$supp, list(testp$jobs), FUN = mean)
colnames(testaggs) <- c("Treatment", "Support")
testaggp <- aggregate(testp$supp, list(testp$price), FUN = mean)
colnames(testaggp) <- c("Treatment", "Support")
testplot <- rbind(testaggc, testaggs, testaggp)
testplot$Treatment[testplot$Treatment=="white collar"] <- "White collar"
testplot$Treatment[testplot$Treatment=="blue collar"] <- "Blue collar"
testplot$Treatment <- factor(testplot$Treatment, levels = testplot$Treatment)
testplot$num <- round(testplot$Support, 2)

ggplot(testplot, aes(y = Support, x = Treatment, fill = c("red", "red", "blue", "blue", "green", "green"))) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_text(aes(label=num), position=position_dodge(width=0.9), vjust=-0.25, size = 8) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_text(size = 32),
        legend.position = "none", axis.text.x = element_text(size = 30, angle = 30, vjust = 0.6),
        axis.text.y = element_text(size = 30)) +
  coord_cartesian(ylim= c(0.4, 0.7))
  

testaggc <- aggregate(testp$supp, list(testp$partisan_dri), FUN = mean)
colnames(testaggc) <- c("Party", "Support")
testaggc$Party <- c("Democrat", "Republican", "Independent")
testaggc$num <- round(testaggc$Support, 2)
testaggc$Party <- factor(testaggc$Party, levels = testaggc$Party)

ggplot(testaggc, aes(y = Support, x = Party, fill = c("red", "blue", "green"))) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_text(aes(label=num), position=position_dodge(width=0.9), vjust=-0.25, size = 8) +
  theme(panel.grid.minor=element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_text(size = 32),
        legend.position = "none", axis.text.x = element_text(size = 30, angle = 30, vjust = 0.6),
        axis.text.y = element_text(size = 30)) +
  coord_cartesian(ylim= c(0.4, 0.7))

#### Data pre-processing ####
results$natsupagg <- (results$nat_sup_1 + results$nat_sup_2) / 2

#### Regression Modeling ####
#### Main regressions
reg1 <- lm(reshoring ~ I(country=="China") + I(jobs=="blue collar") + I(price==1250), data = results)
reg1a <- lm(reshoring ~ I(country=="China") + I(jobs=="blue collar") + I(price==1250) + 
              natsupagg + gender_male + age + 
              income + education + I(partisan_dri == 2) + trump, data = results)
stargazer(reg1, reg1a, omit.stat = c("rsq", "ser", "f"),
          covariate.labels = c("China", "Blue collar", "High price", "National superiority",
                               "Male", "Age", "Income", "Education", "Republican",
                               "Trump voter"))
# Results presented in Table A4 in the Appendix








### Reviewers requested bar plot for dependent variable, Figure A4 in the Appendix.
# Reviewers asked to change two graphs in the appendix. These graphs were frequency charts of the dependent variables, but presented integer data as smooth. Below, we opt for a simple bar chart instead.
# Create a temporary factor variable for plotting purposes
results <- results %>%
  mutate(reshoring_label = factor(reshoring, levels = 1:5, labels = c("Strongly Disagree", "Disagree", "Neutral", "Agree", "Strongly Agree")))

# Calculate frequencies
freq_data <- results %>%
  count(reshoring_label) %>%
  mutate(percentage = n / sum(n) * 100)

# Define custom colors
custom_colors <- c("Strongly Disagree" = "#636363", "Disagree" = "#969696", "Neutral" = "#bdbdbd", "Agree" = "#d9d9d9", "Strongly Agree" = "#f0f0f0")

# Create bar chart with custom colors
ggplot(freq_data, aes(x = reshoring_label, y = percentage, fill = reshoring_label)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = sprintf("%.1f%%", percentage)), vjust = -0.5) +
  scale_x_discrete(name = "Support for Reshoring Policies") +
  scale_y_continuous(name = "Percentage") +
  scale_fill_manual(values = custom_colors, name = "Support Level") + # Custom subdued colors
  theme_minimal() +
  theme(legend.position = "none")
ggsave("/Users/USER/FOLDER/SUBFOLDER/Reshoring/dvdist4_replication.pdf")







```

```{r Summary Statistics}
### We want to make a summary stats table with all the interesting IPE covariates for Appendix figures A1 and A2.

## Pretreatment Ones
summ_df <- results %>%
  select(`Support` = reshoring, 
         `Ideology` = partisan_scale_dr, 
         `Male` = gender_male, 
         `House Reps Term Correct?` = pol_knowledge_2, 
         `Knowledge of COMPETES` = pol_knowledge_3, 
         `Work with UN` = coop_inter_1, 
         `Work with Others` = coop_inter_2, 
         `Defend Human Rights` = coop_inter_3, 
         `Help Developing World` = coop_inter_4, 
         `Solve Conflicts Worldwide` = iso_1, 
         `Americans First` = iso_2, 
         `Out of World Affairs` = iso_3, 
         `Interests of US First` = nat_1, 
         `War for US` = nat_2, 
         `National Attitude` = nat_att_1, 
         `US Superior` = nat_sup_1, 
         `Ashamed of US` = nat_sup_2, 
         `Trust Feds` = gov_trust, 
         `Trust Other Countries` = int_trust, 
         `iPhone Sourcing Correct?` = iphone, 
         `Globalization Good?` = correctglob)


summ_df <- summ_df %>% mutate(Male = case_when(Male == 0 ~ "Female",
                                               Male == 1 ~ "Male"))



datasummary(Support +
              `Ideology` + 
         `House Reps Term Correct?`+ 
         `Knowledge of COMPETES`+ 
         `Work with UN`+ 
         `Work with Others`+ 
         `Defend Human Rights`+ 
         `Help Developing World`+
         `Americans First`+ 
         `Out of World Affairs`+ 
         `Interests of US First`+ 
         `War for US`+ 
         `National Attitude`+ 
         `US Superior`+ 
         `Ashamed of US`+
         `Solve Conflicts Worldwide`+
         `Trust Feds`+ 
         `Trust Other Countries`+ 
         `iPhone Sourcing Correct?`+ 
         `Globalization Good?` ~ as.factor(Male) * (Mean + SD + Histogram), data = summ_df, output = "SummTable3.html")

# If in docx instead.
datasummary(Support +
              `Ideology` + 
         `House Reps Term Correct?`+ 
         `Knowledge of COMPETES`+ 
         `Work with UN`+ 
         `Work with Others`+ 
         `Defend Human Rights`+ 
         `Help Developing World`+
         `Americans First`+ 
         `Out of World Affairs`+ 
         `Interests of US First`+ 
         `War for US`+ 
         `National Attitude`+ 
         `US Superior`+ 
         `Ashamed of US`+
         `Solve Conflicts Worldwide`+
         `Trust Feds`+ 
         `Trust Other Countries`+ 
         `iPhone Sourcing Correct?`+ 
         `Globalization Good?` ~ as.factor(Male) * (Mean + SD + Histogram), data = summ_df, output = "SummTable3.docx")


```

```{r Wave 2 HTE Graphs}
### First, we want to see whether the effect is driven by race; and by manufacturing sector workers. This is for Figure 3 in the main text of the paper.

# Manufacturing

testp <- testp %>% mutate(
    
    manufacturer = case_when(industry == 1 ~ "Non-Manufacturing",
                             industry == 2 ~ "Non-Manufacturing",
                             industry == 3 ~ "Non-Manufacturing",
                             industry == 4 ~ "Non-Manufacturing",
                             industry == 5 ~ "Non-Manufacturing",
                             industry == 6 ~ "Non-Manufacturing",
                             industry == 7 ~ "Non-Manufacturing",
                             industry == 8 ~ "Non-Manufacturing",
                       industry == 9 ~ "Manufacturing",
                       industry == 10 ~ "Non-Manufacturing",
                       industry == 11 ~ "Non-Manufacturing",
                       industry == 12 ~ "Non-Manufacturing",
                       industry == 13 ~ "Non-Manufacturing",
                       industry == 14 ~ "Non-Manufacturing",
                       industry == 15 ~ "Non-Manufacturing",
                       industry == 16 ~ "Non-Manufacturing",
                       industry == 17 ~ "Non-Manufacturing",
                       industry == 18 ~ "Non-Manufacturing",
                       industry == 19 ~ "Non-Manufacturing",
                       industry == 20 ~ "Non-Manufacturing"))

grouped_data <- testp %>%
  select(country, partisan_dri, country_therm, supp, partisan_scale_dr, globalization, price, nostalgia, gender_male, unemployed, education, socio_ingroup, gov_trust, Region, us_citizen, percent, jobs, coop_inter_2, pol_knowledge_2, iphone, manufacturer) %>%
  group_by(manufacturer, country, jobs) %>%
  summarize(mean_support = mean(supp), se_support = sd(supp)/sqrt(n()))

pfg1 <- ggplot(grouped_data, aes(country, mean_support, fill = jobs, group = jobs)) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_errorbar(aes(ymin=mean_support-se_support, ymax=mean_support+se_support), width=.2,
                 position=position_dodge(.9), color = "red4") +
  #geom_point(aes(color = as.factor(jobs)), position = position_dodge(0.01),size=2, shape = 21, fill = "white") +
  scale_fill_manual(values = c("blue4","grey")) +
  #geom_line() +
  labs(
    x = "Country Treatment",
    y = "Support",
    color = "Reshored Job Type",
    title = "HTE by Manufacturing Sector"
   ) +
  theme_bw() +
  theme(text = element_text(size = 10)) +
  labs(fill = "Job Type") +
  #scale_fill_discrete(name = "Reshored Job Type", labels = c("Blue Collar", "White Collar")) +
  #ylim(2.5,4.5) +
  #geom_hline(yintercept=c(3,4), color = "red", size = 0.3) +
  facet_grid(~manufacturer)
pfg1

    
ggsave("/Users/USER/FOLDER/SUBFOLDER/Reshoring/ManufacturingHTE_replication.pdf", 
       width = 7, height = 7)





# Race
testp <- testp %>% mutate(
    
    WhiteNonwhite = case_when(race == 1 ~ "White",
                             race == 2 ~ "Non-White",
                             race == 3 ~ "Non-White",
                             race == 4 ~ "Non-White",
                             race == 5 ~ "Non-White",
                             race == 6 ~ "Non-White",
                             race == 7 ~ "Non-White"))


grouped_data <- testp %>%
  select(country, partisan_dri, country_therm, supp, partisan_scale_dr, globalization, price, nostalgia, gender_male, unemployed, education, socio_ingroup, gov_trust, Region, us_citizen, percent, jobs, coop_inter_2, pol_knowledge_2, WhiteNonwhite) %>%
  group_by(WhiteNonwhite, country, jobs) %>%
  summarize(mean_support = mean(supp), se_support = sd(supp)/sqrt(n()))

fg4 <- ggplot(grouped_data, aes(country, mean_support, fill = jobs, group = jobs)) +
  geom_bar(position = 'dodge', stat = "identity") +
  geom_errorbar(aes(ymin=mean_support-se_support, ymax=mean_support+se_support), width=.2,
                 position=position_dodge(.9), color = "red4") +
  #geom_point(aes(color = as.factor(jobs)), position = position_dodge(0.01),size=2, shape = 21, fill = "white") +
  scale_fill_manual(values = c("blue4","grey")) +
  #geom_line() +
  labs(
    x = "Country Treatment",
    y = "Support",
    color = "Reshored Job Type",
    title = "HTE by Race"
   ) +
  theme_bw() +
  theme(text = element_text(size = 10)) +
  labs(fill = "Job Type") +
  #scale_fill_discrete(name = "Reshored Job Type", labels = c("Blue Collar", "White Collar")) +
  #ylim(2.5,4.5) +
  #geom_hline(yintercept=c(3,4), color = "red", size = 0.3) +
  facet_grid(~WhiteNonwhite, labeller = labeller(unemployed = 
    c("0" = "White",
      "1" = "Non-White"), ncol=2))
fg4

ggsave("/Users/USER/FOLDER/SUBFOLDER/Reshoring/RaceHTE.pdf", 
       width = 7, height = 7)





```

```{r Chips in the News}
### Chips were in the news around the time of our surveys. 
### We will use the "news" dataset from above and then add vertical lines where the surveys were run to make figure A1 in the Appendix.

# Convert to class date
news <- news %>% 
  mutate(DATE = as.Date(DATE, format = "%m/%d/%y"))
news$STEEL


scale_color_manual(name='Regression Model',
                     breaks=c('Linear', 'Quadratic', 'Cubic'),
                     values=c('Cubic'='pink', 'Quadratic'='blue', 'Linear'='purple'))

# Make graph
news %>% 
  ggplot(aes(x=DATE)) +
    geom_line(aes(y = SEMICONDUCTORS), color="#69b3a2") +
    geom_line(aes(y = STEEL), color="yellow4") +
    #ylim(0,1300) +
    scale_x_date(date_labels = "%Y", name = "Year") +
    geom_vline(aes(xintercept = as.Date(DATE[144])), linetype = 4, colour = "red4") +
    geom_vline(aes(xintercept = as.Date(DATE[153])), linetype = 4, colour = "blue1") +
    geom_text(aes(x=news$DATE[140], y= 4500, label = "Study 1"), colour = "red4", vjust = 1, size = 3, angle = 90) +
    geom_text(aes(x=news$DATE[149], y= 4500, label = "Study 2"), colour = "blue1", vjust = 1, size = 3, angle = 90) +
    #scale_x_continuous(name="Years") +
    scale_y_continuous(name="Number of Articles", limits = c(0,5000)) +
    ggtitle("News Articles Over Time Involving Term 'Semiconductor' Versus 'Steel'") +
    theme_classic() +
    geom_text(aes(x=news$DATE[15], y= 500, label = "Semiconductors"), colour = "#69b3a2", vjust = 1, size = 4) +
    geom_text(aes(x=news$DATE[10], y= 1500, label = "Steel"), colour = "yellow4", vjust = 1, size = 4) +
    theme(text=element_text(size=14))

   ggsave("/Users/USER/FOLDER/SUBFOLDER/Reshoring/SemiconductorNews_replication.pdf", 
           width = 9, height = 5)



   
### For the other appendix Google trends graphs, Figures A2 and A3 in the Appendix respectively.
goochips <- read.csv("trends_chips.csv")
goochips$week <- as.Date(goochips$week, format = "%m/%d/%y")
ggplot(goochips, aes(y = num, x = week)) +
  geom_line() +
  geom_vline(xintercept=goochips$week[148],linetype=2, color="red") +
  geom_text(x = goochips$week[155], label="Study 1", colour="red", angle=90, y = 90) +
  geom_vline(xintercept=goochips$week[192],linetype=2, color="blue") +
  geom_text(x = goochips$week[199], label="Study 2", colour="blue", angle=90, y = 90) +
  labs(x = "Week", y = "Salience") +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
        axis.title.x = element_text(size = 32), 
        axis.title.y = element_text(size = 32),
        axis.text.x = element_text(size = 30),
        axis.text.y = element_text(size = 30))

goosteel <- read.csv("trends_steel.csv")
goosteel$week <- as.Date(goosteel$week, format = "%m/%d/%y")
ggplot(goosteel, aes(y = num, x = week)) +
  geom_line() +
  geom_vline(xintercept=goosteel$week[148],linetype=2, color="red") +
  geom_text(x = goosteel$week[155], label="Study 1", colour="red", angle=90, y = 90) +
  geom_vline(xintercept=goosteel$week[192],linetype=2, color="blue") +
  geom_text(x = goosteel$week[199], label="Study 2", colour="blue", angle=90, y = 90) +
  labs(x = "Week", y = "Salience") +
  theme(panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
        axis.title.x = element_text(size = 32), 
        axis.title.y = element_text(size = 32),
        axis.text.x = element_text(size = 30),
        axis.text.y = element_text(size = 30))
```






